Multilayer Graphene as an Endoreversible Otto Engine

We examine the performance of a finite-time, endoreversible Otto heat engine with a working medium of monolayer or multilayered graphene subjected to an external magnetic field. As the energy spectrum of multilayer graphene under an external magnetic field depends strongly on the number of layers, so too does its thermodynamic behavior. We show that this leads to a simple relationship between the engine efficiency and the number of layers of graphene in the working medium. Furthermore, we find that the efficiency at maximum power for bilayer and trilayer working mediums can exceed that of a classical endoreversible Otto cycle. Conversely, a working medium of monolayer graphene displays identical efficiency at maximum power to a classical working medium. These results demonstrate that layered graphene can be a useful material for the construction of efficient thermal machines for diverse quantum device applications.


Introduction
There is rapidly growing interest in the development of quantum technologies, devices that take advantage of the unique properties of quantum systems to enhance their performance. This, in turn, has led to an increased focus on the field of quantum thermodynamics [1][2][3][4]. Within the broad spectrum of topics that fall under the umbrella of quantum thermodynamics, significant attention is paid to the study of quantum heat engines-devices that extend the principles of classical heat engines to include working mediums made up of quantum systems [5][6][7][8][9].
When considering possible systems to use as the working medium of a quantum heat engine, graphene stands out as an intriguing candidate. Graphene's optical, electronic, and mechanical properties have been extensively studied in recent years [10][11][12][13][14][15][16]. Furthermore, graphene has been shown to be useful in a wide range of applications spanning everything from high-accuracy sensing [17][18][19][20][21], to use as an adsorbent for environmental cleanup [22]. This motivates the question as to whether graphene can also be applied to develop highly efficient or powerful nanoscale thermal machines. Monolayer graphene is particularly notable in this context, since its low energy excitations behave as relativistic massless fermions [23]. As such, the study of quantum heat engines with graphene as a working medium can give insight into the role of relativistic quantum features in thermodynamic behavior [24,25].
When graphene layers are stacked to form a multilayer system, the behavior of the system can vary depending on the stacking order. For example, the relativistic character of the charge carriers is only preserved for some layer numbers and stacking arrangements, leading to a unique energy spectrum when in the presence of a perpendicular magnetic field [26]. As the applied magnetic field produces quantized Landau levels in the energy spectrum, the thermodynamic response of these materials presents strong dependence on the stacking type and number of graphene layers [27]. In the context of thermodynamic devices, it was previously shown that the performance of a quasistatic Otto engine with a working medium of twisted bilayer graphene subjected to a magnetic field achieves the highest efficiency when the twist angle of the bilayer corresponds to the magic angle of 0.96 degrees [28]. This latter result demonstrates that the number and configuration of the graphene sheets plays a significant role in the engine performance. Significant attention has also been given to graphene-based engines in the context of continuous, thermoelectric machines [29][30][31][32]. Notably, graphene has also been used in the construction of an experimental nanoscale cyclic heat engine [33].
Efficiency, defined as the net work per cycle divided by the heat absorbed into the working medium from the hot thermal reservoir, is by far the most prominent metric of engine performance. Of similar physical relevance is the power output, defined as the ratio of the net work per cycle to the total cycle time To maximize engine efficiency, the strokes of the cycle must be carried out quasistatically. However, truly quasistatic strokes require an infinite cycle duration, thus leading to vanishing power output. Practically useful metrics of heat engine performance must, therefore, account for cycles implemented in finite time. Endoreversible thermodynamics [34][35][36] provides a framework for introducing finite-time behavior by assuming that, while the working medium remains in a state of local equilibrium at all times during the cycle, the heating and cooling strokes occur quickly enough that the working medium never fully thermalizes with the hot and cold reservoirs. A prominent performance characteristic within endoreversible thermodynamics is the efficiency at maximum power (EMP), which corresponds to maximizing the power output with respect to the external control parameter and then determining the efficiency at that maximum power output. Endoreversible cycles have also been studied in the context of quantum heat engines. We draw particular attention to Ref. [37], where it was shown that the EMP of an endoreversible Otto cycle with a quantum harmonic oscillator as the working medium exceeds the Curzorn-Albhorn (CA) efficiency [34], the EMP achieved by the Otto cycle with a classical working medium.
The heat engine performance analysis carried out in Ref. [28] assumes quasistatic operation in which the bilayer graphene working medium always reaches thermodynamic equilibrium at reservoir temperatures. In this work, we expand the study of cyclic graphene heat engines into two notable directions. First, we examine how the number of graphene layers in the working medium impacts the engine performance, determining analytical expressions for the efficiency and power output valid for monolayer and multilayer systems. We also demonstrate how the number of layers impacts the parameter regimes where the cycle functions as an engine versus a refrigerator. Second, we go beyond the quasistatic assumption and analyze the finite time performance of the engine using the framework of endoreversible thermodynamics. This allows us to examine additional performance metrics, such as the EMP, which we compare to a standard benchmark for EMP, the CA efficiency.
In Section 2, we provide relevant background, including the analytical results of the energy spectrum for monolayer, bilayer, and trilayer graphene. In Section 3, we determine a closed form for the partition function and examine the equilibrium thermodynamic behavior of multilayer graphene. In Section 4, we introduce the endoreversible Otto cycle for multilayer graphene before presenting the results for the engine efficiency, power output, EMP, and parameter regimes where the cycle functions as an engine or refrigerator in Section 5. Finally, in Section 6, we provide some discussion of our results in the context of experimental implementations before concluding in Section 7.

Model
In multilayer graphene structures, the graphene sheets are placed on top of each other in different stacking configurations and are connected through weak van der Waals forces. The stacking configurations are determined by the orientation of the two triangular sublattices that make up the primary honeycomb lattice of a single sheet. For two stacked sheets, three possible orientations, A, B, and C, are possible, each corresponding to displacing one of the sublattice atoms along the edge of the honeycomb with respect to the neighboring sheet [38]. Subject to a perpendicular external magnetic field, these systems can be analyzed using a π-orbital continuum model. Such an analysis is described extensively in ref. [38]. In our analysis we will focus on two particular stacking configurations. For bilayer graphene, we consider Bernal stacking, also known as AB stacking. For the case of trilayer graphene, we consider the rhombohedral configuration, also known as ABC stacking [39]. Significantly, an exact analytical result for the energy spectrum as a function of the external magnetic field can be found for these two cases as we show below.

Monolayer Graphene
In monolayer graphene, the application of a perpendicular magnetic field results in unevenly spaced Landau levels with an energy spectrum proportional to the root of the level quantum number n [40], where B is the magnitude of the magnetic field, e is the electron charge,h is Planck's constant, and v f is the Fermi velocity (∼ 10 6 m/s). Such an energy spectrum is characteristic of ultra-relativistic massless particles with the Fermi velocity playing the role of the speed of light. In Figure 1a, we show the Landau levels for monolayer graphene, where the positive energy curves correspond to electrons and the negative ones to holes [38]. Each one of these Landau levels is four-fold degenerate, including the zero energy state, where the factor of four arises from spin degeneracy and non-equivalent BZ points K and K , known as valley degeneracy.

Bilayer Graphene: AB Stacking
For a bilayer system, the most stable coupling corresponds to Bernal or AB stacking. This consists of displacing the A sublattice atoms of the upper layer so that they lie on top of the B sublattice atoms of the lower layer. Contrary to the monolayer case, the bilayer system has a quadratic dispersion relation, which gives rise to an interesting phenomenon. While the Dirac equation still models the dynamics of the low energy states, the quadratic dispersion relation indicates that the described charge carriers have mass. In this case, under a perpendicular external magnetic field, the Landau level spectrum for the bilayer graphene presents a linear dependence with the magnetic field, as shown in Figure 1b, which takes the form [41], E n = ±hω c n(n − 1), n = 0, 1, 2, . . . , where ω c ≡ eB/m * corresponds to the cyclotron frequency. This effective mass, m * , is related to the Fermi velocity and the interlayer interaction parameter, t ⊥ , by m * = t ⊥ /2v 2 f . This corresponds to a numerical value of m * ∼ (0.039 ± 0.002) m e , where m e is the electron rest mass. Note that Equation (4) has two zero energy levels, corresponding to n = 0 and n = 1. Each of these is eight-fold degenerate, where the factor of eight accounts for the spin, valley, and layer degeneracy.

Trilayer Graphene: ABC Stacking
Rhombohedral ABC-stacked trilayer graphene is one of the most stable crystals commonly obtained in experimental procedures. In the pristine case, it can be considered a zero-gap semiconductor material. When external electric-field potentials are applied to it, it behaves as a semiconductor [42]. The energy spectrum with an external magnetic field has the form [38], Figure 1c shows the Landau levels obtained from the previous Equation (5). In this case, the zero-energy state is twelve-fold degenerate, accounting for spin, valley, and threelayer degeneracy, while the other energy states remain four-fold degenerate, just as in the bilayer and monolayer case.

Partition Function and Equilibrium Thermodynamics
Comparing the energy spectra presented in the previous section, we see a pattern emerging in how the energy scales with the magnetic field. Each energy is proportional to B N /2 , where N is the number of layers. To illustrate this behavior, we plot the first ten positive and negative energy states as a function of the external field for monolayer, bilayer, and trilayer graphene in Figure 1. We also see the energy spectra follow a common structure in regard to the quantum number n, which takes the form, Therefore, we can compactly write the energy spectrum for the multilayer system in the form [38], where is a constant that depends on the number of layers and the stacking structure of the system.
Note that the energy spectra in Equations (3)-(5) include both positive and negative energy solutions corresponding to electrons and holes, respectively. It is important to highlight that our approximation considers the energy spectra close to the Dirac points of the Brillouin zone, as additional energy solutions require extended approaches that we have not yet fully solved. We focus on the electron excitations only, similar to the established analysis for the case of a twisted bilayer graphene system [28]. This regime can be achieved by gating the monolayer graphene [43] and through moderate electron doping for the case of bilayer [44] and trilayer graphene samples [45]. This makes our restriction to the positive energy solutions experimentally available for the three graphene systems we study.
To accurately determine the partition function, we need to carefully consider the degeneracy of the energy levels, especially for the zero-energy state. The compact form of counting these degenerate states in the partition function is given by, The partition function for the energy spectrum given in Equation (7) does not have a simple closed-form solution, except for the case of one layer. However, if we assume that the number of layers N is not very large compared with the number of states n, then the energy spectrum can be approximated as For large n, we can approximate the partition function sum as an integral of the form, Noting that, we obtain a simple analytical form for the partition function, The equilibrium thermodynamic properties of the multilayer graphene system can be determined from the partition function as follows, where F is the free energy, S is the entropy, U is the internal energy, C B is the heat capacity, and M is the magnetization. Most relevant for the performance analysis of an Otto cycle are the internal energy, which is used to determine the work output of the cycle, and the entropy, which must be held constant during the isentropic strokes of the cycle. In Figures 2 and 3, we compare the internal energy (U ) and entropy (S), respectively, for monolayer, bilayer, and trilayer graphene. To ensure that our analytical approximation for the partition function is valid, we also plot the internal energy and entropy determined from numerical calculations of the partition function sum up to 50,000 terms. We find that the analytical approximation matches very well with the numerical calculation, only differing slightly in behavior as T → 0. Examining Figure 2, we see that, for all three cases, the internal energy increases linearly with T in the high temperature regime. From Figure 3 we see that, consistent with the third law of thermodynamics, the entropy approaches a constant value as T → 0.

The Endoreversible Otto Cycle
The Otto cycle consists of four strokes, illustrated graphically in Figure 4 using an entropy (S)-magnetic field (B) diagram. The first stroke (A → B) is an isentropic compression in which the external field is varied from B 1 to B 2 , while the working medium is isolated from the thermal reservoirs. During this stroke, an amount of work, W comp , must be supplied to compress the working medium. The second stroke (B → C) is an isochoric heating stroke in which the working medium draws an amount of heat, Q in , from the hot reservoir while the external field is held constant. The third stroke (C → D) is an isentropic expansion where the working medium is again disconnected from the thermal reservoirs and the external field is varied from B 2 back to B 1 . During this stroke, an amount of work, W exp , is extracted from the expansion of the working medium. The final stroke (D → A) is an isochoric cooling stroke in which the working medium expels and amount of heat, Q out to the cold reservoir while the external field is held constant. Note that the work parameter (B) plays the role of an inverse volume, increasing during the compression stroke (A → B) and decreasing during the expansion stroke (C → D). . Entropy (S) versus external field (B) diagram for the Otto Cycle. Note that the system is only in contact with the thermal reservoirs during the isochoric (vertical) strokes. Note that in the endoreversible framework, the working medium does not fully thermalize to the temperatures T h and T l of the hot and cold reservoirs at points C and A, respectively. Here, Q in is the amount of heat drawn from the hot reservoir during the heating stroke (B → C) and Q out is the amount of heat expelled to the cold reservoir during the cooling stroke (D → A). Similarly, W comp is the amount of work supplied to the working medium during the compression stroke (A → B), while W exp is the amount of work extracted from the working medium during the expansion stroke (C → D).
Characteristic of the framework of endoreversibility, we will assume the working medium remains in a state of local equilibrium throughout the cycle, but, due to finite-time thermalization strokes, never achieves global equilibrium with the reservoirs. The thermodynamic equation of state for the internal energy of the working medium at each corner of the cycle must thus be expressed in terms of the corresponding temperature, T A , T B , T C , or T D , and the external field strength, B 1 or B 2 . Ultimately, we want to determine expressions for the engine performance figures of merit solely in terms of the experimentally controllable parameters, namely the temperatures of the thermal reservoirs, T l and T h , the magnetic field strengths B 1 and B 2 , and the durations of the heating and cooling strokes, τ h and τ l . In order to do so, we must model the thermal conduction during the isochoric strokes and apply the constraint that the entropy remains constant during the isentropic strokes. For this endoreversible analysis, we will follow the procedure established in ref. [37].
During the isentropic compression stroke (A → B), the working medium is decoupled from the thermal reservoirs. As such, all change in the working medium's internal energy can be associated with work, During the isochoric heating stroke (B → C), the external field is held constant. Thus, the difference in internal energy can be associated entirely with heat, As mentioned above, unlike in the quasistatic case, T C = T h since the working medium does not fully thermalize with the hot reservoir. As the heating stroke is now carried out in finite time, we must determine how the temperature of the working medium changes during the duration of the stroke. The temperatures T B and T C , corresponding to the temperature of the working medium at the beginning and ending of the heating stroke, respectively, must satisfy the conditions, where τ h is the duration of the heating stroke. Consistent with the assumptions of endoreversibility, we model thermal conduction in the working medium using Fourier's law.
In this case, the temperature change from T B to T C can be found by applying Newton's law of cooling, where α h is a constant that depends on the thermal conductivity and heat capacity of the working medium. Solving Equation (17) yields, Just as in the compression stroke, the work extracted during the isentropic expansion stroke (C → D) is found from, During the isochoric cooling stroke (D → A), the heat exchanged with the cold reservoir is given by, where, in analogy to the heating stroke, T A and T D satisfy the conditions, We again apply Fourier's law and Newton's law of cooling to model the temperature change during the stroke, which after solving yields, With expressions for the work conducted and heat exchanged during each stroke of the cycle, we can now determine the net work output per cycle, where, by convention, the work completed by the system is treated as negative.
By definition, the entropy remains constant during the isentropic strokes. We can use this fact to obtain a relationship between the initial and final temperatures and magnetic field strengths during the isentropic strokes. Using dS(T, B) = 0, we obtain the following first-order differential equation, Taking the partial derivatives of the entropy found from Equation (13), we arrive at, Solving Equation (26) for the compression stroke we find, Similarly, solving Equation (26) for the expansion stroke we have, This relationship between the temperature, external field, and number of layers can be seen graphically in Figure 5, where we have plotted curves of constant entropy as a function of the temperature and external field for monolayer, bilayer, and trilayer graphene. (c) Figure 5. Isentropic curves as a function of the temperature and external field for (a) monolayer, (b) bilayer, and (c) trilayer systems. Darker shading indicates lower entropy. Here, we have dimensionless parameters with k B = 1.

Efficiency
We are now in a position where we can determine all of our characterizations of engine performance in terms of experimentally controllable parameters. First, combining Equation (1) with Equations (14), (15), (19), (27), and (28), we arrive at a simple expression for the engine efficiency, We note that this expression is strikingly similar to the classical expression for the Otto efficiency, with the layer number N playing the role of the ratio of heat capacities.

Power
Similarly, by combining Equation (2) with (14), (19), (27), and (28) we arrive at an analytical expression for the power output, where we have defined, Here we have defined κ ≡ B 1 /B 2 as the cycle compression ratio. Note that τ h and τ l are the durations of the heating and cooling strokes, respectively, while γ is a multiplicative factor that implicitly incorporates the duration of the isentropic strokes [37]. Examining Equation (30), we see that the power will vanish under the condition that Σ = Λ. This occurs under three conditions. The first is that κ N /2 → T l /T h . We see from Equation (29) that this corresponds to the limit of Carnot efficiency, under which we would expect the power to vanish. The second and third conditions are when exp (α h τ h ) → 1 and exp (α l τ l ) → 1, respectively. These conditions correspond to the limits of instantaneous thermalization strokes or vanishing thermal conductivity, both of which would prevent heat transfer and thus result in zero power. We also note that the power vanishes in the quasistatic limit of τ l + τ h → ∞. From Equations (18) and (23), we see that this limit yields T 3 = T h and T 1 = T l , which in turn maximizes the efficiency. This is a demonstration of the well-established trade-off between efficiency and power.
The efficiency and power are plotted as a function of the compression ratio, κ, in Figure 6, for the monolayer, bilayer, and trilayer systems. We see that the monolayer system has the lowest efficiency but highest power output, while the opposite is true for the trilayer system. The efficiency and power of the bilayer system falls between the monolayer and trilayer results.

Efficiency at Maximum Power
Due to the inherent trade-off between efficiency and power mentioned above, efficiency alone does not provide the most practically useful metric of engine performance. Instead, this role is played by the efficiency at maximum power. In this case, the EMP is found by maximizing the power output with respect to the compression ratio, κ, and then determining the efficiency corresponding to this value of κ. For a classical Otto cycle, the EMP is given by the Curzon-Ahlborn efficiency [34,46], The CA efficiency is especially notable as, like the Carnot efficiency, it depends only on the temperatures of the reservoirs and not only the characteristics of the working medium. However, unlike the Carnot efficiency, the CA efficiency does not provide an upper bound on the efficiency at maximum power, and can be exceeded under certain circumstances.
Due to the complicated expression for power in Equation (30), we maximize the power numerically. The EMP as a function of the ratio of bath temperatures is shown in Figure 7. We see that for the monolayer case, the EMP is identical to the Curzon-Ahlborn efficiency. This result can be confirmed analytically by taking the derivative of Equation (30) and confirming that it vanishes for N = 1 and κ = T l /T h .
For the bilayer and trilayer systems, however, we see that the EMP exceeds the CA efficiency. The EMP is largest for the bilayer system, decreasing slightly in the trilayer case. This trend continues, with the EMP of larger layer numbers converging back towards the CA efficiency. However, it is important to note that if we increase the number of graphene layers significantly beyond the trilayer case, the assumptions made in determining the closed form of the partition function begin to break down. It has been previously shown for both classical and quantum working mediums that, within the regime of linear response, EMP is bounded by the CA efficiency [47][48][49][50][51][52]. To achieve higher EMP requires going beyond the linear regime or by breaking time-reversal symmetry [53][54][55][56][57][58][59][60][61]. For a cyclic engine, the regime of linear response occurs near the equilibrium limit T l ≈ T h .
To probe the behavior of the EMP for a multilayer graphene working medium in and around the linear response regime, we define T l ≡ T and T h ≡ T. In Figure 8, we plot the EMP for the monolayer, bilayer, and trilayer working mediums in comparison with the CA efficiency for = 1.1, 2, and 10. As expected, for the monolayer system, we see that at all examined values of , the EMP is identical to the CA efficiency. For the bilayer and trilayer working mediums, we see that at = 1.1, close to the equilibrium limit, the EMP is identical to CA, consistent with the results in the works mentioned above. As we move away from the equilibrium limit by increasing we see that, at low bath temperatures, the EMP exceeds CA but that as the temperature increases, the EMP converges back to the CA efficiency. As observed in Figure 7, in the low temperature regime, the bilayer working medium EMP exceeds the CA efficiency by a greater amount than the trilayer working medium. However, as temperature increases, the bilayer EMP converges to CA faster than the trilayer EMP.
From these results, we see that two conditions must be met for the EMP to exceed CA. First, the difference in bath temperatures must be sufficiently far from the equilibrium limit of T l ≈ T h . Second, the cycle must be operating in the low temperature, quantum regime, which, for the multilayer graphene working medium, is determined by the condition θ N B 2 /k B T l 1. This second condition is consistent with the results shown in ref. [37] for harmonic working mediums. Here, T l ≡ T, T h ≡ T l . Parameters are α l = α h = τ l = τ h = 1 and B 2 = 2.

Engine vs. Refrigerator
For any arbitrary choice of parameters, it is not guaranteed that the Otto cycle will function as an engine. In general, there are four possible types of thermal machines, corresponding to all possible combinations of directions of heat and work flow consistent with the first and second laws of thermodynamics. An engine corresponds to positive work output, along with heat flow from the hot bath into the working medium, and from the working medium into the cold bath. A refrigerator corresponds to negative work output, along with heat flow from the cold bath into the working medium and from the working medium into the hot bath. A heater corresponds to negative work output, along with heat flow from the working medium into both baths. Finally, an accelerator corresponds to negative work output along with heat flow from the hot bath into the working medium and from the working medium into the cold bath.
By examining the signs of Equations (14), (15), (19), and (20) across the parameter space, we can determine the regions where the cycle will function as each type of thermal machine. Note, both the heater and accelerator are fundamentally nonequilibrium devices, and thus, we would not expect to find regions of parameter space corresponding to these devices under the assumption of endoreversible behavior.
In Figure 9, we show the regions where the cycle functions as either an engine or a refrigerator as a function of the hot bath temperature and compression ratio. We see that over the same region of parameter space, the layer number has a significant impact on the apportionment between engine and refrigerator. In the monolayer case, we see that the majority of the examined region corresponds to the engine regime, while in the trilayer case, the opposite is true, with a larger portion of the explored space corresponding to the refrigerator regime. The origins of this behavior can be understood from the plot of the power in Figure 6. As the layer number decreases, we see that the reduced power output leads to a smaller region where the work is positive and thus a reduced engine regime. 3.

Discussion
Experimental implementation of the multilayer graphene engine requires precise control of the number of layers as well as a tunable magnetic field of sufficient strength to induce Landau quantization. Fine control over multilayer structures has been demonstrated by folding monolayer graphene nanoribbons [62] and precise multilayer thickness measurements can be accomplished with electron spectroscopy [63]. Strong external magnetic fields can be generated either by direct application [64] or through strain-induced psuedo-magnetic fields [65,66].
Graphene is also a prominent Dirac material, a class of systems whose electronic structure gives rise to charge carriers that behave as relativistic fermions. In particular, the energy spectrum of monolayer graphene can be mapped to that of the relativistic Dirac oscillator [67]. We note that our results for the monolayer working medium are consistent with previous work examining the performance of an endoreversible Otto engine with a relativistic oscillator as the working medium [25].

Conclusions
We have examined the performance of an endoreversible Otto cycle with a working medium of a multilayer graphene system. We have found that all examined performance metrics, including the efficiency, power, EMP, and parameter regions under which the cycle functions as an engine or refrigerator all depend significantly on the number of layers. Most notably, we have found that the EMP for bilayer and trilayer graphene working mediums exceeds the Curzon-Ahlborn efficiency. Conversely, we have demonstrated that the EMP of a monolayer graphene working medium is identical to the CA efficiency. We have also found two conditions necessary for the graphene EMP to exceed the CA efficiency, namely the cycle must be sufficiently far from the equilibrium limit and must be operating in the low temperature regime corresponding to θ N B 2 /k B T l 1. These results may be expanded in future work by utilizing Fermi-Dirac statistics to account for the role of holes in the thermodynamic behavior within the temperature range where the Dirac model for graphene remains valid.